Monte Carlo studies of supersymmetric matrix quantum mechanics 
with sixteen supercharges at finite temperature 
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We present the first Monte Carlo results for supersymmetric matrix quantum mechanics with 
sixteen supercharges at finite temperature. The recently proposed non-lattice simulation enables us 
to include the effects of fermionic matrices in a transparent and reliable manner. The internal energy 
nicely interpolates the weak coupling behavior obtained by the high temperature expansion, and 
the strong coupling behavior predicted from the dual black hole geometry. The Polyakov line takes 
large values even at low temperature suggesting the absence of a phase transition in sharp contrast 
to the bosonic case. These results provide highly non-trivial evidences for the gauge/gravity duality. 
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Introduction. — In the last decade, we have witnessed 
the increasing importance of large- N gauge theories in 
theoretical particle physics. For instance, the holo- 
graphic principle, which was inspired originally by the 
Bekenstein-Hawking formula for the black hole entropy, 
has now been given a concrete manifestation as a conjec- 
tured duality between the strongly coupled large- iV gauge 
theory and the weakly coupled supergravity. The best 
understood example is the AdS/CFT correspondence [H, 
but there are numerous extensions to non-conformal field 
theories as well. In particular, large-A^ gauge theories 
in low dimensions have been studied intensively at fi- 
nite temperature, which revealed intriguing connections 
to the black hole thermodynamics 0, y, 0, S @ • 

Large-iV gauge theories in low dimensions also play 
an important role in formulating superstring/M theories 
non-perturbatively based on the idea of matrix models, 
which was successful in the case of non-critical strings. 
For instance, it is conjectured that critical string/M 
theories can be formulated in terms of matrix models, 
which can be formally obtained by dimensionally reduc- 
ing U(iV) super Yang-Mills theory in ten dimensions to 
D = 0, 1, 2 dimensions. In particular, the D = 1 case [7] 
corresponds to the M Theory 0], which is a hypothet- 
ical eleven-dimensional theory proposed to understand 
the dualities among all the known superstring theories in 
ten dimensions. 

In order to confirm these conjectures or to make use 
of them, it is clearly important to study large-iV gauge 
theories from first principles. Monte Carlo simulation is 
expected to be very useful for that purpose. In particu- 
lar, totally reduced models [9| (th e D = case) have been 



studied in refs. [13, EH, H El Q ■ The complex Pfaffian 



which appears from integration over the fermionic matri- 
ces, causes a technical obstacle in numerical simulation, 
which may be overcome by a new method proposed in 



ref. [141 ]. In fact, the phase of the Pfaffian is speculated 
to induce the spontaneous breaking of SO(10) symmetry 
down to SO (4), a scenario for the dynamical generation 
of four-dimensional space-time 15| suggested from the 
Gaussian expansion method [lfij . 

In the D — 1 case, some sort of "discretization" is 
needed in order to put the theory on computer. Given 
the well-known problems with the conventional lattice 
discretization, three of the authors (M.H., J.N. and 
S.T.) have proposed a non-lattice simulation method [Tt| . 
which is useful for studying supersymmetric quantum 
mechanics. The crucial point was that the gauge symme- 
try is almost trivial in Id, and therefore, we can fix the 
gauge completely. This allows us to introduce a Fourier 
mode cutoff A without violating the gauge symmetry. 
In the bosonic case the new method reproduced the lat- 
tice results in the continuum limit. In the SUSY case 
with four supercharges, it reproduced the results of the 
high temperature expansion in the continuum. The same 
model has been studied in ref. [l8| by the lattice approach 
using a simple lattice action and a more complicated lat- 
tice action preserving half of SUSY. 

In this work we apply the non-lattice simulation 
method to the most interesting case with sixteen super- 
charges. We have reduced the computational effort con- 
siderably by introducing pseudo-fermions based on the 
idea of the Rational Hybrid Monte Carlo (RHMC) algo- 
rithm [19( in the way described in ref. [l8| ■ 

As discussed in ref. [13], our action is nothing but the 
gauge-fixed action in the continuum except for having 
a Fourier mode cutoff. Supersymmetry, which is mildly 
broken by the cutoff, is shown to be restored much faster 
than the continuum limit is achieved. In fact, the con- 
tinuum limit is also approached faster than one would 
naively expect from the number of degrees of freedom. 
This is understandable from the fact that the modes 
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above the cutoff are naturally suppressed by the kinetic 
term. A further (albeit technical) advantage of our for- 
mulation is that we can implement the Fourier acceler- 
ation, which eliminates the critical slowing down com- 
pletely , without extra cost since we are dealing with 
Fourier modes directly. We consider that all the theoreti- 
cal and technical merits of the present approach compen- 
sate the superficial increase in the computational effort 
by the factor of 0(A) compared to the lattice approach 
with the same number of degrees of freedom. 

Simulation techniques. — The model can be obtained 
formally by dimensionally reducing lOd super Yang-Mills 
theory to Id. The action is given by 

S = ±JdttrS.±(D t X i f-±[X i ,X j f 

+ ^l> a Dti>a-^i>ahi)al3[Xi,i>i3i\ , (1) 

where D t — dt — i[A(t), ■ ] represents the covariant 
derivative with the gauge field A(t) being an N x N Her- 
mitian matrix. It can be viewed as a one-dimensional 
U(iV) gauge theory with adjoint matters. The bosonic 
matrices Xj(t) (i = 1, • • • ,9) come from spatial compo- 
nents of the lOd gauge field, while the fermionic matri- 
ces ip a (t) ( a — 1) ' ' ' > 16) come from a Major ana- Weyl 
spinor in lOd. The 16 x 16 matrices 7 M in (TT]) act on 
spinor indices and satisfies the Euclidean Clifford alge- 
bra {ji , jj } — 2Sij . We impose periodic and anti-periodic 
boundary conditions on the bosons and fermions, respec- 
tively. The extent j3 in the Euclidean time direction then 
corresponds to the inverse temperature [3 = 1/T. The 
parameter g in (TT]) can always be scaled out by an appro- 
priate rescaling of the matrices and the time coordinate 
t. We take g = without loss of generality. 

We take the static diagonal gauge A(t) = 
■gdiag(ai, • • • , ajv), where a a can be chosen to satisfy the 
constraint max a (a a ) — min a (a a ) < 2ir by using the large 
gauge transformation with a non-zero winding number. 
We have to add to the action a term 
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a<b 



(2) 



which appears from the Faddeev-Popov procedure, and 
the integration measure for a a is taken to be uniform. 
We make a Fourier expansion 

A A 

xf(t) = J2 x%^ nt ; <® = E &- rt , (3) 
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where A = A — 1/2. The indices n and r take integer 
and half-integer values, respectively, corresponding to the 
imposed boundary conditions. Introducing a shorthand 
notation 
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we can write the action ([!]) as S = Sh + Sf , where 
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The fermionic action St may be written in the form 
St = \M.Aar;Bf)si> J tr^%i where we have expanded tp ar = 

^ar^ A m terms of U(AT) generators t A . Integrat- 
ing out the fermionic variables, one obtains the Pfaffian 
PfAl, which is complex in general. However, we ob- 
serve that it is actually real positive with high accuracy 
in the temperature regime studied in the present work. 
Hence we can replace it by |PfA4| = det(I? 1 / 4 ), where 
V = M^M. 

The trick of the RHMC algorithm is to use the rational 
approximation x" 1 ^ ~ 6o + Ylk=i x +b k > wn i cn nas suffi- 
ciently small relative error within a certain range required 
by the system to be simulated. (The real positive param- 
eters ak and bk can be obtained by a code [2l| based on 
the Remez algorithm.) Then the Pfaffian is replaced by 
\PfM\ = J dFdF*e- SpF , where 



<Spf = boF* F ■ 
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using the auxiliary complex variables F, which is called 
the pseudo-fermions. 

We apply the usual HMC algorithm to the whole sys- 
tem as described in ref. [13], except that now we intro- 
duce the momentum variables conjugate to the pseudo- 
fermions F as well as the bosonic matrices Xi and the 
gauge variables a a . When we solve the auxiliary classi- 
cal Hamiltonian dynamics, it is important to apply the 
Fourier acceleration [2(| to the pseudo-fermions F and 
the bosonic matrices JQ. The main part of the compu- 
tation comes from solving a linear system (T> + bk)x 
I (k = 1, • • • , Q). We solve the system for the smallest 
afe using the conjugate gradient method, which reduces 
the problem to the iterative multiplications of A4 to a 
pseudo-fermion field, each of which requires 0(A 2 A^ 3 ) 
arithmetic operations if implemented carefully. The solu- 
tion for larger fe^'s can be obtained as by-products using 
the idea of the multi-mass Krylov solver 22j . This avoids 
the factor of Q increase of the computational effort. 

Infrared instability. — Since the integration domain 
for the bosonic matrices is non-compact, the conver- 
gence of the partition function is not obvious. In par- 
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ticular, there exists a potential danger in the flat di- 
rection corresponding to commuting matrices. Such an 
issue has been addressed in the totally reduced model 
ToL 11 , [HI, [23[ . In the present D = 1 case, let us expand 
the cutoff theory ([5]) around the commuting background 
XiQ = diag(xa, • • • -iXm) and consider the effective ac- 
tion for the moduli parameters xi a and a a . When both 
T and all of \xi a — Xi b \ are large, the fluctuations become 
very massive, and the one-loop approximation is justified. 
We can easily obtain 

(7) 

where we have defined a ab = &a — ct b and x ab = 
\J {%ia — Xib) 2 ■ In eq. ([7]) we have omitted terms inde- 
pendent of x ab , which actually vanish in the A — > oo 
limit. When x ab <C T and a ab <C 2n, the n = term 
dominates and yields a logarithmic attractive potential 
Wi_ioop ^ L a <fc 41 °g{( a ab) 2 + {(3x ab ) 2 } among a a and 
among Xi a . This agrees with the well-known result in 
the bosonic IKKT model, which describes the high tem- 
perature limit of the present model. In fact one obtains 
x ab ~ T 1 / 4 according to the high temperature expansion 
(HTE) [H]. On the other hand, when T <C x ab <C 2ttAT, 
the denominator and the numerator in eq. (|7|) cancel each 
other almost completely. This implies the existence of an 
instability. 

As T is lowered, the instability region approaches the 
peak at x ab ~ T 1 / 4 representing the high temperature 
behavior. However, since we have a sum over all the 
pairs of indices (a, b), a tiny difference between the de- 
nominator and the numerator is enhanced by the factor 
of N 2 . Then it follows that the lower edge of the insta- 
bility region gets multiplied by N. We found empirically 
that the instability can be avoided by taking N > 
which is consistent with the above considerations. 

At low temperature, one also has to worry about the 
finite A effects. In the case of energy, they are negligible 
for A > |r, whereas for the other observables studied in 
this paper, we need twice as large A. 

Results. — Fig. Q] shows the average of the "extent of 
space" R 2 — J dttr (Xi) 2 . The instability mentioned 
above can be probed by the divergence of this quantity. 
We obtain stable results at sufficiently large N. 

In fig. [5] we plot the absolute value of the Polyakov 
line, which is the order parameter for the SSB of the 
U(l) symmetry. Unlike in the bosonic case [2^, |26| |. 
where the SSB occurs around T ~ 1, the Polyakov line 
is not small even at low T. This implies that there is no 
phase transition in the SUSY case, as predicted by the 
gauge/gravity correspondence d, @. We find that the 
Polyakov line can be fitted nicely to 
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FIG. 1: The "extent of space" is plotted against temperature. 
The dashed line represents the result obtained by HTE up to 
the next leading order for N = 8 [24|. 
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FIG. 2: The absolute value of the Polyakov line is plotted 
against temperature. The dashed line represents the result 
obtained by HTE up to the next leading order for N = 8 [24| . 
The dotted line represents a fit to eq. Q with a = 0.15 and 
b = 0.072. 



In fig. [3] we plot the internal energy defined by E = 
where T is the free energy of the system. In 
practice, we calculate it using a formula, which follows 
from a simple scaling argument [l8| . In our case it reads 
E = -3T [(5 b ) - §{(2A + 1)N 2 - 1}] . Our results in- 
terpolate nicely the weak coupling behavior — calculated 
by the HTE up to the next leading order 24j — and the 
strong coupling behavior E/N 2 — 7.4 • T 2 8 predicted by 
the gau ge /g ravity duality [2| from the dual black-hole ge- 
ometry [27j. The power-law behavior sets in at T ~ 0.5, 



a characteristic behavior in a deconfined theory. 



which is reasonable since the effective coupling constant 
is given by A e ft = 1/T 3 in our convention. 

In ref. [J] the Gaussian expansion method was ap- 
plied to the present model, and the energy obtained 
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at the leading order was fitted nicely to the power law 
E/N 2 = 3.2 ■ T 2 - 7 within 0.25 < T < 1. Their results 
are in reasonable agreement with our data at T ~ 1, but 
disagree at lower temperature. 
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FIG. 3: The energy is plotted against temperature. The 
dashed line represents the result obtained by HTE up to the 
next leading order for TV = 8 [24j]. The solid line represents 
the energy predicted at small T by the gauge/gravity duality. 
The upper left panel zooms up the region, where the power- 
law behavior sets in. 

Summary. — In this paper we have presented the first 
Monte Carlo results for the maximally supersymmetric 
matrix quantum mechanics, which is expected to play a 
very important role in string/M theories. The recently 
proposed non-lattice simulation together with the RHMC 
algorithm enabled us to study the low temperature be- 
havior, which was not accessible by the high tempera- 
ture expansion. As we lower the temperature, we ob- 
served the infrared instability, which was found to be 
eliminated, however, by increasing N. We gave a natural 
explanation to this phenomenon based on the one-loop 
effective action. Our data for the internal energy asymp- 
tote nicely to the result obtained from the dual geom- 
etry, which we consider as a highly nontrivial evidence 
for the gauge/gravity duality in the non-conformal case. 
In particular, our results suggest that the maximally su- 
persymmetric matrix quantum mechanics exactly repro- 
duces not only the power but also the coefficient of the 
power-law behavior obtained from the dual black-hole 
geometry. 
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